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Simulating leaky integrate and fire 
neuron with integers 

A. K. Vidybida* 


Abstract 

The leaky integrate and fire (LIF) neuron represents standard neu¬ 
ronal model used for numerical simulations. The leakage is imple¬ 
mented in the model as exponential decay of trans-membrane voltage 
towards its resting value. This makes inevitable the usage of machine 
floating point numbers in the course of simulation. It is known that 
machine floating point arithmetic is subjected to small inaccuracies, 
which prevent from exact comparison of floating point quantities. In 
particular, it is incorrect to decide whether two separate in time states 
of a simulated system composed of LIF neurons are exactly identical. 
However, decision of this type is necessary, e.g. to figure periodic dy¬ 
namical regimes in a reverberating network. Here we offer a simulation 
paradigm of a LIF neuron, in which neuronal states are described by 
whole numbers. Within this paradigm, the LIF neuron behaves ex¬ 
actly the same way as does the standard floating point simulated LIF, 
although exact comparison of states becomes correctly defined. 

Keywords. Leaky integrate and fire neuron; Floating point calculations; 
Simulations 


1 Introduction 


Leaky integrate and fire (LIF) neuron, iKnightI (119721) : ISteinI fll967(l is a 
universally-accepted “standard” neuronal model in theoretical neuroscience. 
Usa ge of LIF as a m od el allowed obta i ning numerous t heoretical results, see 
e.g. Bnrkitt (200^ b); Lanskv ( 1984 1: Tuckwell ( 1988 1. 

When input impulses are absent, the membrane potential of the LIF 
model decays exponentially. In computer simulations, this makes inevitable 
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the usage of machine floating point arithmetic. However, the floa ting point 


calcuk t ions can be c onsiderably inaccurate in some cases, see e.g. iGoldberg 
(1991); Haves f 2nn,*lh . These cases do not so frequently occur. Therefore, 
if one studies neuronal firing statistics, when a simulation is repeated many 
times with slightly different parameters, the errors due to floating point pit- 
falls can be negligibly small. 

The problem arises if one needs to check if two neuronal states obtained 
in the course of simulation are identical. This kind of check is necessary, 
e.g. to fi g ure p eriodic dynamical regimes of a reverberating network, see 
Vidvbidai (1201111 for an example. In order to check if a neural net is in 


a periodic regime, one needs to check whether two states the net passes 
through at two distinct moments of time are exactly the same. If states of 
individual neurons the net is composed of are described by floating point 
variables, then this expects checking that two floating point numbers are 
exactly equal — the operation, which is not permitted in programming due 
to small, but inevitable inaccuracy adhered to the machine floating point 
arithmetic. Instead, the option is to check w hether the dista nce between two 


numbers is less than a ’’machine epsilon”, iHighaml (120021) . This does not 


help to hnd periodic regimes reliably. Indeed, in this comparison paradigm, 
two different states may happen to be treated as identical and this may 
bring about a fake periodic regime. Taking into account that rounding mode 
depends on operating system and on machine architecture, it will be difficult 
to treat in a consistent manner sustained behavior of a dynamical system 
prone to instability (as do neuronal nets) if simulation is made in floating 
point numbers. 

This problem is avoided in Vidvbid^ ( 2011 1. where the binding neuron 
model was used as neuronal model. The binding neuron states are naturally 
described by integers, and machine arithmetic of integers is perfectly accu¬ 
rate. This admits a routine check whether two numbers are exactly equal, 
and finally, whether two states of a neuronal net are identical, provided those 
states are expressed in terms of whole numbers, exclusively. It would be nice 
to have the same possibility when a LIF model is used for simulations. 

In traditional computer simulations of LIF model driven with input stream 
of stimulating impulses, a dynamical system is simulated, which has states 
described with machine floating point numbers. We denote this dynamical 
system as fpLIF. The purpose of this paper is to offer an approximation, 
which replaces the fpLIF with intLIF — a dynamical system whose states 
are expressed in terms of integers. In what follows, we describe the approx¬ 
imation itself and check its quality by running both fpLIF and intLIF with 
the same stream of input impulses randomly distributed in time. It appears 
that it is possible to choose the approximation parameters in such a way. 
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that dynamics of both fpLIF and intLIF are exactly the same, if considered 
in terms of spiking moments. 


2 Methods 

2.1 The fpLIF model 

The simplest LIF model is considered. Namely, the nenronal state at moment 
t is described as membrane voltage, V{t). The resting state is defined as 
V{t) = 0. Any inpnt impnlse advances membrane potential by h, where 
h> 0. Between two consecutive input impulses, V{t) decays exponentially: 

V{t + s)= 

where s > 0, r — is the membrane relaxation time. Suppose that the 
threshold value for the LIF neuron is Vq. The neuron fires a spike every time 
when V{t) > Vq, and V{t) is reset to zero after that. The set of possible 
values of V(t) is the following interval: 

V{t)e[0;Vo[ (2.1) 

The above mentioned properties of the LIF model can be routinely coded 
with V(t) and h declared as floating point quantities. In this case, possible 
values of V(t) will be those machine floating point numbers, which fall into 
the interval fl2.ll) . And this gives the fpLIF dynamical system. 

2.2 The intLIF model 

2.2.1 Pure decay dynamics 

In numerical simulation of a dynamical system, the time is advanced in dis¬ 
crete steps, having duration of a small fixecQ time-step, dt. This gives an 
approximation of the continuous time t with discrete moments: 

T = {0, dt, 2 dt,3 dt,...} (2.2) 

Due to this fact, the membrane voltage, V{t), also changes in a discrete 
manner from step to step. As a result, in a single run of the pure decay 
dynamics (without input stimulation) the V{t) will pass through only some 
discrete values, missing the intermediate ones. Those discrete values can be 

^See Discussion for a remark about adjustable dt. 


3 




chosen as an approximation of the continuous interval fl2.ip . It is clear that 
the set of the discrete values mentioned depends on the initial value of V. 
To be concrete, let us chose the following approximating set: 

Vq, = {a Vo, Vq, Vq, ... }, where a = (2.3) 

This induces the following representation of the interval fl2.ip : 

[0;K)[={0}U [«K);K,[U[«Vo;«K)[U... (2.4) 

Any value V(t) G [0;Vo[ falls for some n into interval [a’^'^^Vo; a^Vo[ and 
we chose its left end as an approximation for that V{t). The error of this 
approximation is less than 


AVh = (1 - «)V)) (2.5) 

Now, if neuron has membrane potential V{t) = V E [a^^^Vo', cWo[ and 
one intends to describe by values from ¥„ its consequent dynamics due to 
pure decay, the following should be done. The V(t) value should be replaced 
with Vd = a^~^^Vo, and the subsequent decay dynamics simply propels V{t) 
through values aVd, a^Vd, a^Vd, ... from Vq. The state of neuron at each 
time step can be labeled with integer n by the following way: The state 
V {t) = a'^^^Vo obtains label n, where n > 0. Now, the decay dynamics can 
be expressed in terms of n. Namely, if at some discrete moment of time f G T 
a state is labeled as n, then the state at the next moment, t + dt, is labeled 
as n + 1. Thus, having in mind the correspondence: n -H- the 

pure decay dynamics becomes as simple as adding unity to the state label at 
each time step. 

For computer simulation purpose, it should be noticed, that the set Vq 
is finite. If the state label n is declared in a program as int, then Vq 
has INT_MAX+1 elements, where INT_MAX is the largest value of int-type 
variable which can be represented in the operating system. Thus, with state 
labels of type int, instead of fl2.3p . one has to chose as Vq the following set 
of INT_MAX+1 elements: Vq = {a Vq, Vq, ..., q/I hT_MAX where the 

value 0 is added to describe the resting state, which is attained just after 
firing. Consequently, (12.41) should be replaced with the following: 

[0; Vo[= [aVo; Vo[U[aVo; «Vo[U • • • U [0; cr^^T.MAX ( 2 . 6 ) 

In a 64-bits OS, INT_MAX = 2147483647. This imposes a limit on possible 
duration of the pure decay evolution represented with int type labels. If the 
time step dt = 0.01 ms, then INT_MAX-(if ~ 350 minutes. Thus, description 
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of neuronal state by int type label fails, if a LIF neuron starts, e.g. with 
V (t) close to Vo and does not receive excitatory stimulation for longer than 
350 minutes of real time. 

In real networks, stimulus-free period of a neuron involved in a useful task 
cannot be so long. Indirectly, a very crude upper bound for possible duration 
of stimulus-free period can be estimated from duration o f supp ression of 
activity observed in some brain networks, lOssandon et al.l (1201111 . Such a 
suppression can last up to 1000 ms. Based on this value, let us expect that 
our neuron receives excitatory input impulse^ of amplitude h with mean 
rate higher than r ~ 0.1 Hz. Before the hrst input impulse, neuron is in 
the state ’’empty” with V = 0. After the hrst input impulse the state 
label 77-1 ~ logQ,(h/Vo) and this value is small if compared with INT_MAX. 
For example, for dt = 0.01 ms, r = 20 ms, H/Vq = 0.01 one obtains rii ~ 10^. 
The label ni gets increment 1 after each time step dt until the next input 
impulse comes. The mean waiting time for the next stimulation, if expressed 
in the dt units, is l/(r ■ dt) = 10®. This means that at the moment of next 
stimulation, the state label n 2 ~ ni-|-10® <C 10^. After the next input impulse, 
the state label < rii. This suggests that states with labels n > 10”^ will be 
observed quite rare with no chances to attain n value close to the INT_MAX 
provided a neuron receives at least a moderate stimulation. 


2.2.2 Input impulse 

Suppose that input impulse at moment t advances by h the membrane voltage 
V{t). If V{t) G Vq, then V{t) + h ^ Ya does not hold in most cases. One 
needs to approximate V{t) + hhj a. suitable value from Vq,, as described in n. 
12.2.11 and then to proceed with pure decay dynamics expressed in integers, 
as described in n. 12.2.11 until the next input impulse. 

Preliminary tests of this scheme where performed with LIF neuron stim¬ 
ulated with Poisson stream of input impulses. The standard floating point 
LIF (fpLIF) simulation and integer LIF (intLIF) simulation, as described 
here, were performed with the same input streams. The threshold manner 
of triggering gives chances that bring moments of both models will be the 
same, if expressed in whole number of dt. Computer simulations performed 
show that the bring moments of both models are indeed the same for some 
initial period of time after which they become diberent. This is because 
the approximation error, when approximating LIF state (voltage) by a value 
from Ya, builds up too fast with each input impulse. In order to decrease 
the error, one needs more precise approximation of the continuum set [0; Vo[, 

^See the next section for exact treatment of input impulses. 
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than that given by the discrete set Vq,. The reqnired approximation can be 
achieved by introdncing second order bins into the set Vq,. 

2.2.3 Second order bins 

Let us choose an integer > 1 and divide each bin from fl2.6p to N equal 
second order bin^. This gives representation of n-th first order bin: 

N-l 

K+Vo; a^Vo[= IJ [a"+Vo + z • c„; + (* + l)c„[ 

i=0 

where c„ denotes the size of second order bin within the n-th hrst order bin 

Cn = {a^Vo - a"+Vo)/iV (2.7) 

Now, if V{t) G -|- i • c„; -|- (f -|- l)c„[, then we chose Vd = 

a" +Vo + * • c„ as its approximation. This results in the new set Ya,N of 
possible values for V{t): 

INT_MAX-i N-i 

V„,v = {0}U IJ lJ{a"+Vo + zc4 

n=0 2=0 

Any point, except of 0, in the set Ya,N is labelled with two integers, {n, z}, 
where 0 < n < INT_MAX, 0 < z < iV. Any label {n, i} corresponds to the 
membrane voltage 14 ,i from Ya,N- Vn,i = a^+^Vo + tCn, or: 

= + ( 2 . 8 ) 

From the last, it is clear that pure decay evolution for time step dt transforms 
voltage 14 ,i into 14 + 1 ,*, which means for the integer labels: 

{n, z} {n + l,i} (2.9) 

Now, if initial state of neuron is given as floating point number V G [0; 14 [, 
then one needs to find its approximation 14,i G Va,N, as described above. The 
precision of this approximation is: 

|y-14,*| < AV^ (2.10) 

^We do not consider the last bin from (EH), because we do not expect the state V will 
ever fall into it, except of just after firing, when V = 0^ exactly. 
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In the approximation of V with I4i j, the value of n should satisfy the 
following relation: <V< a"‘Vo, which gives: 

n = - llog„(Vr)| - 1 (2.11) 

where [x] denotes the integer part of x. Now, with n found, we determine the 
i value from the following relation: + icn <V < + (i + l)c„, 

which gives: 

t= [(V-a^+%)/cn] ( 2 . 12 ) 

Now evolution of LIF state can be expressed in integer numbers as fol¬ 
lows. For initial value of voltage, 1^(0) < Vq, we hnd its integer representa¬ 
tion/approximation {n,i} in accordance with Eqs. fl2.1ip . 02.121) . The pure 
decay evolution then goes as displayed in 02.9p . In order to describe influ¬ 
ence of input impulse with magnitude h on the state {n, i] we calculate the 
voltage Vn^i in accordance to 02.8p . The voltage after receiving input impulse 
becomes Vnew = + h. If Vnew > ho, then neuron produces an output 

impulse and ends in the state “empty” with V" = 0. Otherwise, the new 
integer state, {n',F}, can be found with 02.111) . 02.12p applied to the Vnew 
instead of V. 

The procedures given above dehne the dynamical system intLIF in which 
a neuronal state is described with two integers, {n, i}, with additional unique 
state “empty” attained just after bring. 


3 Results 

3.1 Coding 

Neuronal state is described by three integers: int n,i; char empty; If 
empty == 1 then the neuron is in its resting state with V = 0. If empty == 0 
then V can be calculated from the n,i in accordance with 02.8p . Equation 
02.8p gives the following C-code for calculating V from known integer state 
■[n,i, empty}: 

V = empty ? 0 : pow(al, n)*V0*(al + (double)i/N*(1 - al)); 

where VO stands for Vq and al stands for a. After bring, empty == 1. 
Equation 02.lip gives the following C-code for calculating n: 

n = - floor(log(VO/V)/log(al)) - 1; 

Similarly for the equations 02.7p . 02.12p . 
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3.2 Testing 

For testing the intLIF simulation paradigm, both intLIF and fpLIF models 
were stimulated with the same random stream of impulses. The stream 
with exponential distribution of inter-spike intervals (ISI) was generated with 
random number generator from the GNU Scientihc Librarjfl. Generators of 
three types were used, each with a number of different seeds, see Table [H 

Table 1: Parameters of simulating algorithm 
random number generators MT19937, knuthran2002, tausllS 

number of different seeds for each 10 
testing duration, real time 1 hour 

max value for N Nmax = 10® 

min value for dt dtmin = 0.001 ms 

initial N 10 

initial dt 0.1, 0.01, 0.001 ms 


Each ISI obtained from the generator as floating point number was ap¬ 
proximated with a value from T by rounding to the nearest integer number 
of time steps dt: 

double ISI = gsl_ran_exponential (r, mu); 
int ilSI = rint(ISI/dt); 

ISI = dt*iISI; 

where gsl_ran_exponential (r, mu) — is the generator, r — is the pointer 
to global generator, mu — is the mean inter-spike interval of the exponential 
distribution; mu = 1/A, where A is the intensity of Poisson stream, which 
is obtained from the generator and used as input stimulating stream. This 
way obtained ISIs were used to apply input impulses to both the intLIF and 
fpLIF model within a single program. The program starts with initial values 
for N, dt given in the Tabled] With each input impulse, it is checked whether 
both fpLIF and intLIF react the same way, namely, either both hre, or both 
not hre. If the reaction for some input impulse is not the same, a new value 
for N is chosen by multiplying the current value by 10, and simulation starts 
anew from the beginning. If the maximum value for N is reached (see Table 
dj) and fpLIF and intLIF bring moments still do not the same through the 
whole simulation time, next option is to divide dt by 10. A single run of 
the program is considered successful if during 1 hour of real time all input 
impulses produced the same reaction in both intLIF and fpLIF. 

^See http://www.gnu.org/software/gsl/ 



Several sets of physical parameters, which cover a physiologically reason¬ 
able range, were nsed in the testing, see Table [H 


Table 2: Physical parameters nsed in simulation. 


threshold depolarization, Vq 
input impulse height, h 
membrane relaxation time, r 
input stream intensity, A 


20 mV 

0.25, 0.5, 1., 2., 4., 8., 16. mV 

10., 20., 40. ms 

0.4, 0.8, 1.60, 3.20, 6.40 ms"^ 


The testing was made for any combination of parameters from both Table 
1 and 2. Some characteristic numbers of input and output streams are given 

in Table [21 


Table 3: 
testing. 


Characteristics of streams of impulses during 1 hour of real time 


minimal number of input impulses 
maximal number of input impulses 
minimal number of output impulses 
maximal number of output impulses 


~ 1.4-10® 
~ 23.4-10® 
0 

~ 11.7-10® 


As a result of testing, it was found that for any combination of parame¬ 
ters it is possible to ensure that all bring moments of intLIF and fpLIF are 
identical by choosing proper N and dt values. The decisive factor, which 
determines whether all spiking moments of intLIF and fpLIF coincide, is the 
precision of approximating the interval of possible voltages, [0; Vo[, with dis¬ 
crete values from Ya,N as compared to h. This relative error, 6V, can be 
estimated based on 02.51) and 02.1UI) : 6V = (1 — a)Vo/(Nh), which for small 
^ gives: 5V = {dt Vq)/ {t N h). In the testing performed, it appeared that 
having 6V < 2.0 - 10“^^ guarantees that sequences of spiking moments of 
intLIF and fpLIF are identical. For larger values of 6V, differences between 
the two sequences may appear, which are characterized by a few misplaced 
spikes, up to several hundred and more. 

4 Discussion 

In numerical simulation of a dynamical system adaptive algorithms are nor¬ 
mally used, when the time step dt value is increased/decreased during cal¬ 
culations in order to make calculations faster and more precise. This works 
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perfect if it is necessary to calculate the system’s state at some future mo¬ 
ment of time starting from some initial state. To determine periodic regimes 
in a reverberating network, it is instead necessary to calculate the whole tra¬ 
jectory during some interval of time. In this case, the straightforward way is 
to approximate that interval with equidistant discrete points as in fl2.2jl and 
calculate states of the system in those points. Therefore, paradigm of hxed 
time step is used here for a single neuron. 

Description of neuronal state (voltage) with a pair of integers {n,z} does 
not exempt the intLIF model from using machine floating point numbers. 
Indeed, in Eqs. fl2.8jl . fl2.11j) . fl2.12jl . operations with floating point numbers 
are explicitly involved. Nevertheless, the pure decay evolution, as it is de¬ 
scribed in fl2.9p goes without rounding errors. The underlying to fl2.9p values 
from are always the same for the same {n,i}- With obtained input 

impulse, calculation of Vn^i + h involves a rounding error. The error is imme¬ 
diately cleared while calculating new {n,i} by means of 02.111) . 02.12p . This 
allows describing states of LIF neuron in terms of integers in a consistent 
manner. Namely, different state labels {n,i} correspond to different voltages 
from Ya,N and vice versa. 

We used here the simplest possible model for LIF neuron. It seems, that 
technique offered in nn. 12.2.11 12.2.21 above, can be ex tended to be va lid 
for mor e elab o rated LIF mode ls, like those described in iBnrkittI fj2nn6alJbh : 
Lanskvl ( 1984 ): Tuckwell ( 1988 ). 


5 Conclusions 

The intLIF paradigm for numerical simulation of leaky integrate and hre 
neuron is proposed. In this paradigm, neuronal state is described by two 
integers, {n, i}. The membrane voltage of LIF neuron can be calculated from 
{n, f}, if required. The LIF state change due to both leakage and stimulating 
impulses is expressed exclusively in terms of changing integers {n,i}. The 
intLIF paradigm is compared with the standard fpLIF simulation paradigm, 
where membrane voltage is expressed as machine floating point number, by 
stimulating both models with the same random stream of input impulses 
and registering the spiking moments of both models. It is concluded that 
approximation parameters of intLIF can be chosen in such a way that spiking 
moments of both models are exactly the same if expressed as whole number 
of simulation time step, dt. Description of LIF states by integers gives a 
consistent numerical model, suitable for situations where exact comparison 
of states is necessary. 
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